﻿using Accessibility;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Numerics;
using System.Text;
using System.Threading.Tasks;

namespace AstroDisplayUI
{
    class SolarSystemCalc
    {
        //based on code from http://stjarnhimlen.se/comp/tutorial.html by Paul Schlyter pausch@stjarnhimlen.se
        public static (double,double) Calculate(DateTime dobs, double obs_latitude_degn, double obs_longitude_dege)
        {
            //dobs = local date and time of observation
            //obs_latitude_degn = observation latitude in degrees north (-90 to 90, 39 for Washington DC), for azimuth and altitude
            //obs_longitude_dege = observation longitude in degrees east (-180 to 180, 77 for Washington DC), for approximate UTC correction, not too important

            DateTime dobs_utc = dobs - TimeSpan.FromDays(obs_longitude_dege / 360);//approximate UTC date of observation
            //the above conversion does not take into account time zones or equation of time
            //since this is only used to calculate orbit parameters which change slowly
            //there is a DateTime.ToUniversalTime() function but it uses an implicit time zone
            DateTime d0 = new DateTime(1999, 12, 31);
            //day number
            TimeSpan dayNum = dobs_utc - d0;
            double d = dayNum.TotalDays; //[days] since Dec 31 1999 UST (approximate) used for orbital parameters
            const double dtor = Math.PI / 180;

            double ecl = (23.4393 - (3.563e-7 * d)%360) * dtor;//[rad] obliquity of the ecliptic

            OrbitalElements Sun = new OrbitalElements() { N = 0, i = 0, w = (282.9404 + (4.70935e-5 * d)%360)*dtor, a = 1, e = 0.016709 - 1.151e-9 * d, M = (356.0470 + (0.9856002585 * d)%360)*dtor };
            OrbitalElements Moon = new OrbitalElements() { N = (125.1228 - (0.0529538083 * d)%360)*dtor, i = 5.1454*dtor, w = (318.0634 + (0.1643573223 * d)%360)*dtor, a = 60.2666, e = 0.0549, M = (115.3654 + (13.0649929509 * d)%360)*dtor };


            //compute position of the sun
            //eccentric anomaly
            double E = Sun.M + Sun.e * Math.Sin(Sun.M) * (1 + Sun.e * Math.Cos(Sun.M));
            double xv = Math.Cos(E) - Sun.e;
            double yv = Math.Sqrt(1 - (Sun.e * Sun.e)) * Math.Sin(E);
            AstroPosition sungeoecl = new AstroPosition(xv, yv, 0);
            sungeoecl.C1 += Sun.w;
            AstroPosition sungeoeq = AstroPosition.EclRotate(sungeoecl, ecl);

            //convert further to local coordinates
            double obs_latitude = obs_latitude_degn * dtor; //[rad]
            double obs_localtime = dobs.TimeOfDay.TotalHours * (Math.PI / 12); //[rad] 0 to 2*pi, pi = noon
            //the above local time does not take into account time zones or equation of time
            //therefore the result may correspond to almanac entries up to ~30 minutes earlier or later than the time specified in dobs
            double L = Sun.w + Sun.M; //mean longitude
            double gmst0 = L + Math.PI;
            double sidtime = obs_localtime + gmst0;
            double ha = sidtime - sungeoeq.C1;
            AstroPosition sungeoloc = new AstroPosition() { C1 = ha, C2 = sungeoeq.C2, Radius = sungeoeq.Radius };
            AstroPosition sunaz = AstroPosition.HorRotate(sungeoloc, Math.PI / 2 - obs_latitude);
            double azimuth = sunaz.C1 + Math.PI;
            double altitude = sunaz.C2;
            return (azimuth, altitude);
        }

        class OrbitalElements
        {
            internal double N, i, w, a, e, M;
        }


        class AstroPosition
        {
            public double Radius, C1, C2;
            //public enum Pos { Equatorial, Ecliptic, Azimuth };
            //Equatorial: C1=RA, C2=Decl [rad]
            //Ecliptic: C1=Longitude, C2=Latitude [rad]
            //Azimuth: C1=Asimuth, C2=Altitude [rad]
            //public Pos Type;

            public double[] GetRectangular()
            {
                double[] v = new double[3];
                v[0] = Radius * Math.Cos(C1) * Math.Cos(C2);
                v[1] = Radius * Math.Sin(C1) * Math.Cos(C2);
                v[2] = Radius * Math.Sin(C2);
                return v;
            }

            public AstroPosition() { }
            public AstroPosition(double x, double y, double z)
            {
                Radius = Math.Sqrt(x * x + y * y + z * z);
                double pd = Math.Sqrt(x * x + y * y);
                C1 = Math.Atan2(y, x);
                C2 = Math.Atan2(z, pd);
            }

            public static AstroPosition EclRotate(AstroPosition p1, double ecl)
            {
                double[] sr = p1.GetRectangular();
                double x = sr[0];
                double y = sr[1] * Math.Cos(ecl) - sr[2] * Math.Sin(ecl);
                double z = sr[1] * Math.Sin(ecl) + sr[2] * Math.Cos(ecl);
                return new AstroPosition(x, y, z);
            }

            public static AstroPosition EquRotate(AstroPosition p1, double ha)
            {
                double[] sr = p1.GetRectangular();
                double x = sr[0] * Math.Cos(ha) - sr[1] * Math.Sin(ha);
                double y = sr[0] * Math.Sin(ha) + sr[1] * Math.Cos(ha);
                double z = sr[2];
                return new AstroPosition(x, y, z);
            }

            public static AstroPosition HorRotate(AstroPosition p1, double ha)
            {
                double[] sr = p1.GetRectangular();
                double x = sr[0] * Math.Cos(ha) - sr[2] * Math.Sin(ha);
                double y = sr[1];
                double z = sr[0] * Math.Sin(ha) + sr[2] * Math.Cos(ha);
                return new AstroPosition(x, y, z);
            }
        }

    }
}
